clc; clear; 
cd 'data'

lambda_data = csvread('lambda.csv');
rho_data = csvread('rho.csv');
Y_data = csvread('Y.csv');

cd ..

I = 44;

alpha = 1; 
theta = 2.5; sigma = 3;
YEAR = 15;

lambda = reshape(lambda_data(:,YEAR),I,I);    
rho = 1-reshape(rho_data(:,YEAR),I,I);   
rho_autarky =  max(0,1-(1-diag(rho))./diag(lambda));
 
Y = Y_data(:,YEAR);

K = 2;
grid = [0.05,0.95];

for p = 1:K
   s = grid(p);
   DOM = 1./(1-diag(rho)).*(((1-s)./(1-diag(rho))).^((1-sigma)/theta) - 1);
   DOM = (s>=diag(rho)).*DOM;
   r(:,p) = (1-alpha*(sigma-1)./(sigma.*theta)) + sum(alpha*(theta-(sigma-1))./(sigma.*theta).*DOM,2);
   LAMBDA = lambda.*repmat(Y',I,1)./repmat(Y,1,I)./(1-rho).*(((1-s)./(1-rho)).^((1-sigma)/theta) - 1);
   LAMBDA = max(0,abs(1-eye(I,I)).*LAMBDA);
   INT(:,p) = alpha*(theta-(sigma-1))./(sigma.*theta).*sum(LAMBDA,2);
   %----------------------------------------------------------------------------------------------------%
   DOM =  1./(1-rho_autarky).*(((1-s)./(1-rho_autarky)).^((1-sigma)/theta) - 1);
   DOM = (s>=rho_autarky).*DOM;
   ra(:,p) = (1-alpha*(sigma-1)./(sigma.*theta)) + sum(alpha*(theta-(sigma-1))./(sigma.*theta).*DOM,2);
end


change = diag(lambda).^(1./theta).*ra./(r+INT);
 
WW = [linspace(1,I,I)',change,diag(lambda).^(1./theta)]; 
WW = sortrows(WW,4);

lb = WW(:,2)'; ub = WW(:,3)';
eps = 0.3;
xx = []; yy = [];

for s = 1:I
    
   x = [s-eps, s+eps, s+eps, s - eps];
   
   xx = [xx x'];
   
   y = [lb(s) lb(s) ub(s) ub(s)];
   
   yy = [yy y'];   
   
end

patch(xx,yy,'black', 'FaceAlpha', 0.7,'EdgeColor','none');
xticks(linspace(1,I,I));
xticklabels({"LUX","MLT","IRL","HUN","BEL","LTU","SVK","EST","SVN","NLD","CZE","CYP","BGR","DNK","TWN","AUT","HRV","LVA","SWE","FIN","POL","PRT","ROW","DEU","ROU","NOR","CHE","CAN","GRC","KOR","MEX","FRA","ESP","GBR","TUR","ITA","IDN","RUS","JPN","AUS","IND","BRA","USA","CHN"}); 

%=========================================================================%

%% Numbers in the paper
mean(yy(1,:), 'all')
mean(yy(3,:), 'all')